'''
Created on Apr 8, 2010

@author: oabalbin
'''
import sys
import numpy as np
import matplotlib.pyplot as plt
import matplotlib.pylab as pla
import networkx as nx
import signatures.parsers.read_gene_lists as gpm
import signatures.db.gene_annotations as ga
# import SOAP Client
from SOAPpy import WSDL
from optparse import OptionParser
from collections import defaultdict, deque
from signatures.common.classes import cPathway, cInteraction


'''
class cInteraction():
    """
    Class interaction to be used in the KEGG script
    """
    def __init__(self,pathname,id1,id2,names1,names2,type_int,type_rel):   
        #item1.element_id1, gen1, item1.element_id2, gen2,type_interact,type_relation
        self.pathname = pathname
        self.item1id = id1
        self.item2id = id2
        self.item1names = names1
        self.item2names = names2
        self.type_interact = type_int
        self.type_relation = type_rel


'''


'''
def kegg_dictOfrels():
    
    relations = {
                'PPrel': ['activation',inhibition,
                PPrel    activation_indirect effect,
                PPrel    inhibition_indirect effect,
                PPrel    activation_phosphorylation,
                PPrel    inhibition_phosphorylation,
                PPrel    activation_dephosphorylation,
                PPrel    inhibition_dephosphorylation,
                PPrel    glycosylation,
                PPrel    methylation,
                PPrel    dephosphorylation,
                PPrel    phosphorylation,
                PPrel    ubiquination,
                PPrel    inhibition_ubiquination,
                PPrel    state change,
                PPrel    binding/association,
                PPrel    dissociation,
                PPrel    complex,
                PPrel    compound,
                PPrel    indirect effect,
                PPrel    hidden compound,
                PPrel    missing interaction,
                GErel    expression,
                GErel    repression,
                GErel    indirect effect,
                GErel    missing interaction,
                ECrel    compound,
                ECrel    hidden compound,
                PCrel
                maplink
                 }
'''

def convert_genelist2keggid(genelist):
    """
    it converts a ncibi-id gene list to a list of KEGG IDs which can be used 
    later for the rest of the queries.   
    """
    wsdl = 'http://soap.genome.jp/KEGG.wsdl'
    serv = WSDL.Proxy(wsdl)

    genelist_tmp = ['ncbi-geneid:'+gen for gen in genelist]
    querylist = ",".join(genelist_tmp).replace(',',' ')
    results = serv.bconv(querylist) 
    results = results.split('\t')
    KEGG_list = deque([])
    
    for item in results:
        f=item.find('hsa:')
        if f != -1:
            KEGG_list.append(item)
            
    print KEGG_list
    return KEGG_list


class cNets():
    
    def __init__(self):
        self.theNet= nx.DiGraph()
        self.tga = ga.gene_annotation("localhost", "oabalbin", "oscar", "gene_annotation")   
    def fill_network(self,inputfile):
        """
        It creats an inderected graph, by using the adjancency list provided in the input file. 
        The adjancency list is created using the KEGG databases and It has 6 fields
        """
        
        for line in inputfile:
            line = line.strip('\n')
            fields = line.split('\t')
            
            #node1, node2 = fields[2],fields[3]
            #relation_type=defaultdict(list)
            #relation_type[fields[4]] = fields[5]
            #relation_type = fields[4]+'_'+fields[5]
            # simplistic view: if GErel then you can affect expression, not otherwise 
            #print node1, node2, relation_type
            # not using undefined relations or enzymatic relations
            if fields[2] == 'undefined' or fields[3]=='undefined' or fields[4]=='ECrel':
                continue
            else:
                
                kegg_node=[fields[2],fields[3]]
                #print kegg_node
                entrez_nodes = [g.replace('hsa:','') for g in kegg_node]
                #print entrez_nodes 
                hugo_nodes= self.tga.translate_entrez2hugo(entrez_nodes)
                #print hugo_nodes
                node1, node2 = hugo_nodes[0], hugo_nodes[1] 
                
                if fields[4]=='GErel':
                    relation_type=1.0
                else:
                    relation_type=0.0

                if relation_type==0.0:
                    continue
                self.theNet.add_nodes_from([node1, node2])
                self.theNet.add_edge(node1,node2,relation_type)
    
    def get_pathsfrom_Node2Rest(self,source_node):
        connected_nodes = nx.dfs_preorder(self.theNet, source=source_node, reverse_graph=False)
        #connected_nodes = nx.node_connected_component(self.theNet, source_node)
        return connected_nodes 
    
    def get_shortest_pathsfrom_Node2Set(self,source_node,target_nodes):
        """
        It obtains the nodes in the shorthest path between a source node and all the nodes in the target_nodes list
        """
        paths=[]
        mynodes = self.theNet.nodes()
        for tn in target_nodes:
            if tn in mynodes:                
                print nx.shortest_path(self.theNet, source_node, tn)
                paths.append(nx.shortest_path(self.theNet, source_node, tn))
            else:
                print tn
            
        return paths
    
    
    def remove_nodes_fromNet(self, listofNodes):
        """
        Removes nodes from a network in the following way: if A ==> B ==> C, and B is removed then A ==> C
        """
        terminalnodes=0
        internalnodes=0
        for n in listofNodes:
            successors = self.theNet.successors(n)
            #print 'out and in degree of node '+n
            #print self.theNet.out_degree(n)
            #print self.theNet.in_degree(n)
            if not successors:
                #print 'no successors'
                self.theNet.remove_node(n)
                terminalnodes+=1                        
            else:
                internalnodes+=1
                predecessors = self.theNet.predecessors(n)
                for p in predecessors:
                    #print 'predecessors'
                    #print p
                    for s in successors:
                        #print 'successors'
                        #print s 
                        relation_type = self.theNet[p][n]
                        #print relation_type 
                        self.theNet.add_edge(p,s, relation_type)
                
                self.theNet.remove_node(n)
        print 'The number of terminal nodes = '+str(terminalnodes)
        print 'The number of internal nodes = '+str(internalnodes)    


    def write_edge_list(self, outfile):
        """
        write the adjancecy list including the attributes.        
        """
        for node,nbrdict in self.theNet.adjacency_iter():
            for neigh, edge_property in nbrdict.iteritems(): 
                attr = edge_property.keys()[0]
                attrvals = edge_property[attr]
                edge = [node,neigh,attr,attrvals]
                outfile.write(",".join(edge).replace(',','\t')+'\n')
    



#### Initialization
inputfile = open('/home/oabalbin/projects/erg_regulation/2010_04_09_ERG_regulation_KEGG')
inputfile2 = open('/home/oabalbin/projects/erg_regulation/2010_04_09_ERG2genesChipSeq')
inputfile3 = open('/home/oabalbin/projects/erg_regulation/2010_04_09_ERG2Interactome')
##
outputfile = '/home/oabalbin/projects/erg_regulation/2010_04_10_s_ERG_network_draft.txt' #shortestpath
outputfileX = '/home/oabalbin/projects/erg_regulation/2010_04_10_s_ERG_Allnetwork_draft.txt'
outputfile4 = '/home/oabalbin/projects/erg_regulation/2010_04_10_s_ERG_shortestpathnetwork_draft.txt' #shortestpath
outputfile3 = '/home/oabalbin/projects/erg_regulation/2010_04_10_s_ERG_network_draft2.txt'
outputfile5 = '/home/oabalbin/projects/erg_regulation/2010_04_10_s_ERG_shortestpathnetwork_draft2.txt'
outputfile2 = open('/home/oabalbin/projects/erg_regulation/2010_04_10_s_ERG_network_draft_genes.txt','w')
outputfile8 = open('/home/oabalbin/projects/erg_regulation/2010_04_10_s_ERG_network_draft_genes2hsa.txt','w')
outputfile6 = open('/home/oabalbin/projects/erg_regulation/2010_04_10_s_ERG_shortestpathnetwork_draft_genes.txt','w')
###
tga = ga.gene_annotation("localhost", "oabalbin", "oscar", "gene_annotation")
gp =gpm.geneparser()
genelist_query = list(set(gp.list_of_names(open('/home/oabalbin/projects/erg_regulation/all_genes_in_ETS_classifier.txt'))))
gene_inthe_expMatrix = gp.list_of_names(open('/home/oabalbin/projects/erg_regulation/gene_inthe_expMatrix.txt'))
#genelist=tga.translate_hugo2entrez(genelist_query)
#print genelist 
#gene_in_classifier_module = convert_genelist2keggid(genelist)
'''
gene_in_classifier_module = deque(['hsa:2078', 'hsa:7941', 'hsa:247', 'hsa:5226', 'hsa:1298', 'hsa:9766', 'hsa:2690', 'hsa:4217', 'hsa:7358', 'hsa:776', \
                                   'hsa:1627', 'hsa:4646', 'hsa:5192', 'hsa:4950', 'hsa:26751', 'hsa:3065', 'hsa:397', 'hsa:1889', 'hsa:4035', 'hsa:5152', \
                                   'hsa:8618', 'hsa:5754', 'hsa:5873', 'hsa:128', 'hsa:7088', 'hsa:5796', 'hsa:272', 'hsa:7917', 'hsa:3790', 'hsa:10656', \
                                   'hsa:1050', 'hsa:493', 'hsa:23423', 'hsa:5335', 'hsa:7181', 'hsa:290', 'hsa:2171', 'hsa:8764', 'hsa:347', 'hsa:55884', \
                                   'hsa:10647', 'hsa:65108', 'hsa:81544', 'hsa:221443', 'hsa:3109', 'hsa:4015', 'hsa:2678', 'hsa:6904', 'hsa:3020', 'hsa:10397', \
                                   'hsa:8531', 'hsa:1746', 'hsa:9759', 'hsa:4188', 'hsa:5460', 'hsa:1031'])
'''
'''
#genelistChip = list(set(gp.list_of_names(open('/home/oabalbin/projects/erg_regulation/ERG_chipSeq1000'))))
genelistChip = list(set(gp.list_of_names(open('/home/oabalbin/projects/erg_regulation/erg_complex_genecore'))))
print genelistChip
genelistERGChip=tga.translate_hugo2entrez(genelistChip)
print genelistERGChip
print len(genelistERGChip)
# Open the KEEGG api

outputfile7 = open('/home/oabalbin/projects/erg_regulation/all_genes_in_ETS_classifier2.txt','w')
# Get Kegg gene IDs
KEGG_list = convert_genelist2keggid(genelist)
#KEGG_list.append('hsa:2078')
print KEGG_list

for i in KEGG_list:
    line = ['0','0','hsa:2078',i,'PPrel','potential binding/association']
    #line = ['0','0','hsa:2078',i,'GErel','potential expression']
    print line
    outputfile7.write(i+'\n')
    #outputfile7.write(",".join(line).replace(',','\t')+'\n')

sys.exit(0)
'''

##############################

### Network algorithm
mynet = cNets()
mynet.fill_network(inputfile)
print mynet.theNet.number_of_nodes()
print mynet.theNet.number_of_edges()
#
mynet.fill_network(inputfile2)
mynet.fill_network(inputfile3)

print mynet.theNet.number_of_nodes()
print mynet.theNet.number_of_edges()

#ERG='hsa:2078'

ERG='ERG'

# For now: remove self loops from the network
print nx.is_directed_acyclic_graph(mynet.theNet)
print mynet.theNet.number_of_selfloops()
edges_selfloops = mynet.theNet.selfloop_edges(data=True)
mynet.theNet.remove_edges_from(edges_selfloops)
print mynet.theNet.number_of_selfloops()
print nx.is_directed_acyclic_graph(mynet.theNet)
'''
# Get all connected components in the graph.
Gcc=nx.strongly_connected_component_subgraphs(mynet.theNet)
G0=Gcc[0]
ind=0
print len(G0)
for Gi in Gcc[0:]:
    if len(Gi)>0:
        if ERG in Gi.nodes():
            print len(Gi)
            print 'ERG is here' 
        ind+=2
print ind
'''

# get all the nodes that belong to the ERG component, 
# by using the first depth search algorithm
ERG_component = mynet.get_pathsfrom_Node2Rest(ERG)
# Generate a subgraph using only the nodes in the ERG component
ERG_subgraph = nx.subgraph(mynet.theNet, ERG_component, copy=True)

# How many classifier genes are in the network
nodes_inNet = set(mynet.theNet.nodes())
commonNodes = nodes_inNet.intersection(set(genelist_query))
print 'Number of common Nodes in the net'
print len(commonNodes), len(set(genelist_query))
print commonNodes

erg_successors = mynet.theNet.successors(ERG)
print 'number of ERG succesors'
print len(set(erg_successors))
predecessors=[]
successors=[]
'''
for n in commonNodes:
    if n == 'hsa:1050':
        print mynet.theNet.predecessors(n)
        print mynet.theNet.successors(n)
    predecessors.append(mynet.theNet.predecessors(n))
    successors.append(mynet.theNet.successors(n))
'''

for n in erg_successors:
    successors.append(mynet.theNet.successors(n))

    
predecessors = set(sum(predecessors,[]))
successors = set(sum(successors,[]))

print 'number of prdecessor and successors'
print len(predecessors), len(successors)

commonNodes_neighbors = predecessors.union(successors)
print len(commonNodes_neighbors)
#commonNodes_neighbors = commonNodes_neighbors.union(set(commonNodes))
commonNodes_neighbors.add(ERG)
commonNodes_neighbors = commonNodes_neighbors.union(set(erg_successors))
print len(commonNodes_neighbors)

commonNodes_subgraph = nx.subgraph(mynet.theNet, list(commonNodes_neighbors), copy=True)

print 'number of terminal nodes in the final network'
print commonNodes_neighbors.intersection(commonNodes)

print len(commonNodes_subgraph.nodes())
print len(commonNodes_subgraph.edges())


ETS1 = ['ETS1']
HDAC8 = ['HDAC8']
erg2targetsNet = cNets()
erg2targetsNet.theNet = commonNodes_subgraph
genes_in_net = set(erg2targetsNet.theNet.nodes())
genes_without_expression = genes_in_net.difference(set(gene_inthe_expMatrix))
print 'genes with out expression'
print genes_without_expression
print len(genes_without_expression)
erg2targetsNet.remove_nodes_fromNet(list(genes_without_expression))

nx.write_adjlist(erg2targetsNet.theNet,outputfile)
nx.write_adjlist(mynet.theNet,outputfileX)
nx.write_edgelist(erg2targetsNet.theNet,outputfile3,delimiter='\t')

sys.exit(0)

# Translate the genes from kegg ids back to hugo ids 
#ERG_component_entrez = [g.replace('hsa:','') for g in commonNodes_subgraph]
#ERG_component_hugo = tga.translate_entrez2hugo(ERG_component_entrez)
for gen in commonNodes_subgraph.nodes():
    outputfile2.write(gen+'\n')


sys.exit(0)

'''  
nodes_ERGcomp = set(ERG_subgraph.nodes())
common2ERGcomp = nodes_ERGcomp.intersection(commonNodes)
out2ERGcomp = commonNodes.difference(common2ERGcomp)
print len(common2ERGcomp), len(set(commonNodes))
print common2ERGcomp
print len(out2ERGcomp)
print list(out2ERGcomp)

# Calculate the boundary of a component
nodes_inboundary = nx.node_boundary(mynet.theNet, ERG_component, out2ERGcomp)
edges_inboundary = nx.edge_boundary(mynet.theNet, ERG_component, out2ERGcomp)
print len(ERG_component)
print nodes_inboundary
print len(nodes_inboundary)
print edges_inboundary
# Are the classifier genes in the nodes of the strong component
print nx.shortest_path(mynet.theNet, ERG, 'hsa:1050')
'''
# Draw the network, giant component and other components
'''
try:
    from networkx import graphviz_layout
    layout=nx.graphviz_layout
except ImportError:
    print "PyGraphviz not found; drawing with spring layout; will be slow."
    layout=nx.spring_layout
print 'making layout'
pos=layout(mynet.theNet)
print 'Calculating connected components'
Gcc=nx.connected_component_subgraphs(mynet.theNet)
G0=Gcc[0]
print 'Drawing network'
nx.draw(mynet.theNet,pos, with_labels=False,node_size=10)
nx.draw_networkx_edges(G0,pos,with_labels=False,edge_color='r', width=6.0)
# show other connected components
ind=0
for Gi in Gcc[1:]:
    if len(Gi)>1:
        print ind
        ind+=1
        nx.draw_networkx_edges(Gi,pos,with_labels=False,edge_color='r',alpha=0.3,width=5.0)

plt.savefig("/home/oabalbin/projects/erg_regulation/giant_component.png")
#plt.show() # display


'''

#print ERG_subgraph.neighbors('hsa:861')
#print ERG_subgraph.neighbors('hsa:1050')
#print ERG_subgraph.neighbors('hsa:2078')

#ergNet = cNets()
#ergNet.theNet = ERG_subgraph
#paths2targets = mynet.get_shortest_pathsfrom_Node2Set(ERG,gene_in_classifier_module)
#paths2targets = list(set(sum(paths2targets,[])))
#print paths2targets
#erg2targetes = nx.subgraph(mynet.theNet, paths2targets, copy=True)

#erg2targetsNet = cNets()
#erg2targetsNet.theNet = erg2targetes

# write the network in the adjancecy list format
nx.write_adjlist(ERG_subgraph,outputfile)
nx.write_adjlist(mynet.theNet,outputfileX)
#nx.write_adjlist(erg2targetes,outputfile4)
#nx.write_multiline_adjlist(ERG_subgraph, outputfile, delimiter='\t', comments='#')
#nx.write_multiline_adjlist(erg2targetes, outputfile, delimiter='\t', comments='#')
nx.write_edgelist(ERG_subgraph,outputfile3,delimiter='\t')
#nx.write_edgelist(erg2targetes,outputfile5,delimiter='\t')
#ergNet.write_edge_list(open(outputfile3,'w'))
#erg2targetsNet.write_edge_list(open(outputfile3,'w'))

# Translate the genes from kegg ids back to hugo ids 
ERG_component_entrez = [g.replace('hsa:','') for g in ERG_component]
#ERG_component_entrez = [g.replace('hsa:','') for g in erg2targetes]
ERG_component_hugo = tga.translate_entrez2hugo(ERG_component_entrez)
for gen in ERG_component_hugo:
    outputfile2.write(gen+'\n')

print len(set(ERG_component))
#print len(set(erg2targetes))
print mynet.theNet.number_of_nodes()










'''
genelistChip = gp.list_of_names(open('/home/oabalbin/projects/erg_regulation/All_ERG_chip'))
genelistERGChip=tga.translate_hugo2entrez(genelistChip)
print genelistERGChip
print len(genelistERGChip)
# Open the KEEGG api
outputfile7 = open('/home/oabalbin/projects/erg_regulation/2010_04_08_ERG2genesChipSeq','w')
# Get Kegg gene IDs
KEGG_list = convert_genelist2keggid(genelistERGChip)
KEGG_list.append('hsa:2078')
print KEGG_list

for i in KEGG_list:
    line = ['0','0','hsa:2078',i,'GErel','potential expression']
    print line
    outputfile7.write(",".join(line).replace(',','\t')+'\n')

sys.exit(0)
'''



        
        
        
        